Plot conditional effects and grid predictions from 02

Author

Max Lindmark

Published

September 9, 2023

Intro

Plot conditional effects from predictions in script “02-fit-sdms-conditional.qmd”

Load packages & source functions

# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools", "patchwork",
          "viridis", "RColorBrewer", "here", "ggnewscale", "ggh4x") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

pal_temp <- brewer.pal(n = 10, name = "Paired")[c(2, 6)]
pal_oxy <- brewer.pal(n = 10, name = "Paired")[c(8, 4)]

# Set path
home <- here::here()

Read predictions

preds <- read_csv(paste0(home, "/output/data_conditional_02_sdms.csv")) |> 
  filter(oxy > 3 & oxy < 9) |>
  filter(temp > 2 & temp < 14) |> 
  separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage)) |> 
  group_by(group) |> 
  mutate(est_sc = exp(est)/max(exp(est))) # for plotting all together in heatmap
Rows: 15000 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): group
dbl (14): temp, oxy, quarter, temp_sc, temp_sq, oxy_sc, year, year_f, quarte...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 5,700 rows (38%), 9,300 rows remaining

filter: removed 2,418 rows (26%), 6,882 rows remaining

mutate: changed 6,882 values (100%) of 'species' (0 new NA)

        changed 6,882 values (100%) of 'life_stage' (0 new NA)

group_by: one grouping variable (group)

mutate (grouped): new variable 'est_sc' (double) with 3,621 unique values and 0% NA
grid_preds <- read_csv(paste0(home, "/output/data_pred_grids_02_sdms.csv")) |> 
  separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage))
Rows: 5660568 Columns: 32
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): substrate, month_year, ices_rect, group
dbl (28): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: changed 5,660,568 values (100%) of 'species' (0 new NA)

        changed 5,660,568 values (100%) of 'life_stage' (0 new NA)

Calculate weighted quantiles

wq <- grid_preds |> 
  group_by(group, quarter) |>
  summarise("Median_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.5),
            "10_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.1),
            "25_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.25),
            "75_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.75),
            "90_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.9),
            "Env_oxy" = median(oxy),
            
            "Median_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.5),
            "10_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.1),
            "25_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.25),
            "75_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.75),
            "90_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.9),
            "Env_temp" = median(temp)) |> 
  ungroup() |> 
  separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage)) 
group_by: 2 grouping variables (group, quarter)
summarise: now 12 rows and 14 columns, one group variable remaining (group)
ungroup: no grouping variables
mutate: changed 12 values (100%) of 'species' (0 new NA)
        changed 12 values (100%) of 'life_stage' (0 new NA)
wq_annual <- grid_preds |> 
  group_by(group, year, quarter) |>
  # Could use reframe here... 
  summarise("Median_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.5),
            "1st_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.25),
            "3rd_oxy" = hutils::weighted_quantile(v = oxy, w = exp(est), p = 0.75),
            "Env_oxy" = median(oxy),
            
            "Median_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.5),
            "1st_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.25),
            "3rd_temp" = hutils::weighted_quantile(v = temp, w = exp(est), p = 0.75),
            "Env_temp" = median(temp)) |> 
  ungroup() |> 
  separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage)) 
group_by: 3 grouping variables (group, year, quarter)
summarise: now 348 rows and 11 columns, 2 group variables remaining (group, year)
ungroup: no grouping variables
mutate: changed 348 values (100%) of 'species' (0 new NA)
        changed 348 values (100%) of 'life_stage' (0 new NA)

Plot conditional effects

# Oxygen
preds_oxy <- preds |> 
  filter(temp %in% sort(unique(preds$temp))[0.5*length(unique(preds$temp))])
filter (grouped): removed 6,696 rows (97%), 186 rows remaining
wi_range_oxy <- preds_oxy |> 
  left_join(wq |>
              filter(quarter == 4) |>
              dplyr::select(group, `10_oxy`, `25_oxy`, `75_oxy`, `90_oxy`), by = "group")
filter: removed 6 rows (50%), 6 rows remaining
left_join: added 4 columns (10_oxy, 25_oxy, 75_oxy, 90_oxy)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     186
           >                 =====
           > rows total       186
# TODO: make predictions for the exact value??
wi_median_oxy <- preds_oxy |> 
  left_join(wq |>
              filter(quarter == 4) |>
              dplyr::select(group, Median_oxy), by = "group") |> 
  filter(oxy > 0.985*Median_oxy & oxy < 1.015*Median_oxy) |> 
  as.data.frame()
filter: removed 6 rows (50%), 6 rows remaining
left_join: added one column (Median_oxy)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     186
           >                 =====
           > rows total       186
filter (grouped): removed 180 rows (97%), 6 rows remaining
p1 <- ggplot() + 
  ggh4x::facet_grid2(life_stage ~ species, scales = "free", independent = "y") +
  # NOTE: 80% Confidence interval
  
  # 25th and 75th percentile
  geom_ribbon(data = wi_range_oxy |> filter(oxy > `25_oxy` & oxy < `75_oxy`), fill = "black",
              aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
              color = NA, alpha = 0.2) +
  geom_line(data = wi_range_oxy |> filter(oxy > `25_oxy` & oxy < `75_oxy`), aes(oxy, exp(est)),
            color = "black", alpha = 0.5, linewidth = 1.5) + 
  
  # 10th and 90th percentile
  geom_ribbon(data = wi_range_oxy |> filter(oxy > `10_oxy` & oxy < `90_oxy`), fill = "black", 
              aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
              color = NA, alpha = 0.2) +
  geom_line(data = wi_range_oxy |> filter(oxy > `10_oxy` & oxy < `90_oxy`), aes(oxy, exp(est)),
            color = "black", alpha = 0.5, linewidth = 1) + 
  
  # NOTE: 85% Confidence interval
  geom_ribbon(data = preds_oxy,  fill = "black",
              aes(oxy, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)), color = NA, alpha = 0.2) +
  geom_line(data = preds_oxy, aes(oxy, exp(est)), alpha = 0.5, linewidth = 0.5, color = "black") + 
  
  # Median!
  geom_point(data = wi_median_oxy, aes(Median_oxy, exp(est)),
             color = "black", alpha = 0.5, size = 2.5) + 
  
  labs(x = "Oxygen (ml/L)", y = "Biomass density (kg/km<sup>2</sup>)") +  
  theme(legend.position = "bottom",
        legend.spacing.y = unit(0.1, 'cm'),
        legend.key = element_rect(fill = "grey95", inherit.blank = TRUE),
        axis.title.y = ggtext::element_markdown(),
        plot.margin = unit(c(0,0,0,0), "cm")) +
  NULL
filter (grouped): removed 141 rows (76%), 45 rows remaining
filter (grouped): removed 141 rows (76%), 45 rows remaining
filter (grouped): removed 101 rows (54%), 85 rows remaining
filter (grouped): removed 101 rows (54%), 85 rows remaining
# Temperature
preds_temp <- preds |> 
  filter(oxy %in% sort(unique(preds$oxy))[0.5*length(unique(preds$oxy))])
filter (grouped): removed 6,660 rows (97%), 222 rows remaining
wi_range_temp <- preds_temp |> 
  left_join(wq |>
              filter(quarter == 4) |>
              dplyr::select(group, `10_temp`, `25_temp`, `75_temp`, `90_temp`), by = "group")
filter: removed 6 rows (50%), 6 rows remaining
left_join: added 4 columns (10_temp, 25_temp, 75_temp, 90_temp)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     222
           >                 =====
           > rows total       222
# TODO: make predictions for the exact value??
wi_median_temp <- preds_temp |> 
  left_join(wq |>
              filter(quarter == 4) |>
              dplyr::select(group, Median_temp), by = "group") |> 
  filter(temp > 0.985*Median_temp & temp < 1.015*Median_temp) |> 
  as.data.frame()
filter: removed 6 rows (50%), 6 rows remaining
left_join: added one column (Median_temp)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     222
           >                 =====
           > rows total       222
filter (grouped): removed 216 rows (97%), 6 rows remaining
wi_median_temp
       temp    oxy quarter   temp_sc   temp_sq     oxy_sc year year_f quarter_f
1  7.698927 5.7852       1 0.6633954 0.4400935 -0.2013043 1999   1999         1
2  8.024739 5.7852       1 0.7803021 0.6088714 -0.2013043 1999   1999         1
3 10.305422 5.7852       1 1.5986492 2.5556793 -0.2013043 1999   1999         1
4  9.653798 5.7852       1 1.3648358 1.8627766 -0.2013043 1999   1999         1
5  8.350551 5.7852       1 0.8972089 0.8049837 -0.2013043 1999   1999         1
6  8.350551 5.7852       1 0.8972089 0.8049837 -0.2013043 1999   1999         1
  sal_sc depth_sc depth_sq             group  species life_stage        est
1      0        0        0         cod_adult      Cod      Adult  3.4416694
2      0        0        0      cod_juvenile      Cod   Juvenile  1.4838354
3      0        0        0      plaice_adult   Plaice      Adult -1.4838181
4      0        0        0   plaice_juvenile   Plaice   Juvenile -6.4750500
5      0        0        0    flounder_adult Flounder      Adult  4.8314757
6      0        0        0 flounder_juvenile Flounder   Juvenile -0.5487353
     est_se    est_sc Median_temp
1 0.6249370 0.9116380    7.602427
2 0.8294039 0.6392413    8.129976
3 0.8646872 0.7731394   10.387263
4 1.2343494 0.5440757    9.573396
5 0.3224810 0.4976172    8.387172
6 0.5137997 0.4212632    8.333579
p2 <- ggplot() + 
  ggh4x::facet_grid2(life_stage ~ species, scales = "free", independent = "y") +
  # NOTE: 80% Confidence interval
  
  # 25th and 75th percentile
  geom_ribbon(data = wi_range_temp |> filter(temp > `25_temp` & temp < `75_temp`), fill = "black",
              aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)), 
              color = NA, alpha = 0.2) +
  geom_line(data = wi_range_temp |> filter(temp > `25_temp` & temp < `75_temp`), aes(temp, exp(est)),
            color = "black", alpha = 0.5, linewidth = 1.5) + 
   
  # 10th and 90th percentile
  geom_ribbon(data = wi_range_temp |> filter(temp > `10_temp` & temp < `90_temp`), fill = "black", 
              aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
              color = NA, alpha = 0.2) +
  geom_line(data = wi_range_temp |> filter(temp > `10_temp` & temp < `90_temp`), aes(temp, exp(est)),
            color = "black", alpha = 0.5, linewidth = 1) + 
  
  # NOTE: 85% Confidence interval
  geom_ribbon(data = preds_temp,  fill = "black",
              aes(temp, exp(est), ymin = exp(est - est_se*1.282), ymax = exp(est + est_se*1.282)),
              color = NA, alpha = 0.2) +
  geom_line(data = preds_temp, aes(temp, exp(est)), alpha = 0.5, linewidth = 0.5, color = "black") + 
  
  # Median!
  geom_point(data = wi_median_temp, aes(Median_temp, exp(est)),
             color = "black", alpha = 0.5, size = 2.5) + 
  
  labs(x = "Temperature (°C)", y = "Biomass density (kg/km<sup>2</sup>)") +  
  theme(legend.position = "bottom",
        legend.spacing.y = unit(0.1, 'cm'),
        legend.key = element_rect(fill = "grey95", inherit.blank = TRUE),
        axis.title.y = ggtext::element_markdown(),
        plot.margin = unit(c(0,0,0,0), "cm")) +
  NULL
filter (grouped): removed 178 rows (80%), 44 rows remaining
filter (grouped): removed 178 rows (80%), 44 rows remaining
filter (grouped): removed 138 rows (62%), 84 rows remaining
filter (grouped): removed 138 rows (62%), 84 rows remaining
(p1 / p2) + plot_annotation(tag_levels = "A")

ggsave(paste0(home, "/figures/oxy_temp_conditional.pdf"), width = 17, height = 17, units = "cm")

Plot weighted quantiles over time

t <- wq_annual |> 
  pivot_longer(c("Median_oxy", `1st_oxy`, `3rd_oxy`, "Env_oxy", "Median_temp", `1st_temp`, `3rd_temp`, "Env_temp")) |> 
  mutate(var = ifelse(grepl("oxy", name), "Oxygen", "Temperature"),
         type = ifelse(grepl("Env", name), "Environment", "Biomass-weighted")) |> 
  separate(group, c("species", "life_stage"), sep = "_", remove = FALSE) |> 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage),
         group2 = paste(species, life_stage, sep = " ")) |> 
  filter(name %in% c("Median_oxy", "Median_temp", "Env_oxy", "Env_temp"))
pivot_longer: reorganized (Median_oxy, 1st_oxy, 3rd_oxy, Env_oxy, Median_temp, …) into (name, value) [was 348x13, now 2784x7]
mutate: new variable 'var' (character) with 2 unique values and 0% NA
        new variable 'type' (character) with 2 unique values and 0% NA
mutate: changed 2,784 values (100%) of 'species' (0 new NA)
        changed 2,784 values (100%) of 'life_stage' (0 new NA)
        new variable 'group2' (character) with 6 unique values and 0% NA
filter: removed 1,392 rows (50%), 1,392 rows remaining
t_env <- t |>
  filter(type == "Environment" & group == "cod_adult") |> 
  mutate(group2 = "Environment")
filter: removed 1,276 rows (92%), 116 rows remaining
mutate: changed 116 values (100%) of 'group2' (0 new NA)
t2 <- bind_rows(t_env, t |> filter(!type == "Environment"))
filter: removed 696 rows (50%), 696 rows remaining
# Reorder to have Environment at the end
t2$group2 <- factor(t2$group2, levels = rev(unique(t2$group2)))

pal <- brewer.pal(name = "Dark2", n = 8)[3:8]

t2 |> 
  mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter4")) |> 
  ggplot(aes(year, value, color = group2, linetype = group2)) +
  geom_line(alpha = 0.8) +
  #facet_grid(quarter ~ var) +
  ggh4x::facet_grid2(var ~ quarter2, scales = "free") +
  labs(y = "Environmental variable", x = "Year", linetype = "", color = "") +
  guides(color = guide_legend(nrow = 3,override.aes = list(linetype = c(rep(1, 6), 2), 
                                                           size = 0.5)),
         linetype = "none") +
  scale_x_continuous(breaks = c(seq(1993, 2021, by = 5))) +
  scale_linetype_manual(values = c(rep(1, 6), 2)) +
  scale_color_manual(values = c(rev(pal), "grey40")) +
  theme(legend.position = c(0.24, 0.36),
        legend.spacing.y = unit(-2, 'cm'),
        legend.spacing.x = unit(0, 'cm'),
        legend.background = element_rect(fill = NA))
mutate: new variable 'quarter2' (character) with 2 unique values and 0% NA

ggsave(paste0(home, "/figures/wm_temp_oxy.pdf"), width = 17, height = 14, units = "cm")

Extra plot for Gotje presentation

Correlation and slope plot

#  Correlation plot
dcor <- t |> filter(!type == "Environment") |>
  mutate(id = paste(name, year, quarter, sep = "_")) |>
  dplyr::select(group2, year, var, value, quarter, id)
filter: removed 696 rows (50%), 696 rows remaining
mutate: new variable 'id' (character) with 116 unique values and 0% NA
t_env2 <- t_env |>
  mutate(name = ifelse(name == "Env_oxy", "Median_oxy", "Median_temp")) |>
 mutate(id = paste(name, year, quarter, sep = "_")) |>
  rename(env_value = value) |> 
  dplyr::select(id, env_value)
mutate: changed 116 values (100%) of 'name' (0 new NA)
mutate: new variable 'id' (character) with 116 unique values and 0% NA
rename: renamed one variable (env_value)
dcor2 <- dcor |> left_join(t_env2, by = "id")
left_join: added one column (env_value)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     696
           >                 =====
           > rows total       696
cors <- plyr::ddply(dcor2, c("group2", "var", "quarter"),
                    summarise, cor = round(cor(env_value, value), 2)) |> 
  arrange(var)
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
summarise: now one row and one column, ungrouped
hist(cors$cor)

# Plot correlations between environment and weighted
cor_pal <- brewer.pal(n = 8, name = "Dark2")[6:7]
  
dcor2 |> 
  mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter 4")) |> 
  ggplot(aes(env_value, value, color = quarter2)) +
  ggh4x::facet_grid2(var ~ group2, scales = "free") +
  geom_abline(color = "grey30", linetype = 2, linewidth = 0.5) +
  geom_point(alpha = 0.7, size = 0.5) +
  geom_text(data = cors |> filter(quarter == 1),
            aes(label = paste("r=", cor, sep = "")), x = -Inf, y = Inf, hjust = -0.1, vjust = 1.5,
            inherit.aes = FALSE, fontface = "italic", color = cor_pal[2], size = 2.5) +
  geom_text(data = cors |> filter(quarter == 4),
            aes(label = paste("r=", cor, sep = "")), x = -Inf, y = Inf, hjust = -0.1, vjust = 2.7,
            inherit.aes = FALSE, fontface = "italic", color = cor_pal[1], size = 2.5) +
  labs(y = "Biomass-weighted", x = "Environment", color = "") +
  stat_smooth(method = "lm", se = FALSE, size = 0.5) +
  scale_color_manual(values = cor_pal) +
  scale_x_continuous(breaks = c(4:8)) +
  theme_sleek(base_size = 10) + 
  theme(aspect.ratio = 1,
        legend.key.size = unit(0.2, 'cm'),
        legend.text = element_text(size = 6),
        legend.position = c(0.11, 0),
        legend.spacing.y = unit(-1, 'cm'))
mutate: new variable 'quarter2' (character) with 2 unique values and 0% NA
filter: removed 12 rows (50%), 12 rows remaining
filter: removed 12 rows (50%), 12 rows remaining
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
`geom_smooth()` using formula = 'y ~ x'

ggsave(paste0(home, "/figures/supp/wm_cors.pdf"), width = 17, height = 7, units = "cm")
`geom_smooth()` using formula = 'y ~ x'
# Plot the slopes over time for group and quarter, compare to the environment slope. Boxplots and lines??
s_slopes <- t |> filter(!type == "Environment") |>
  mutate(id = paste(group, var, quarter, sep = "_")) |> 
  ungroup() |> 
  dplyr::select(year, id, value)
filter: removed 696 rows (50%), 696 rows remaining
mutate: new variable 'id' (character) with 24 unique values and 0% NA
ungroup: no grouping variables
t_slopes <- t_env |>
  mutate(id = paste("Environment", group2, var, quarter, sep = "_")) |> 
  dplyr::select(year, id, value)
mutate: new variable 'id' (character) with 4 unique values and 0% NA
slopes <- bind_rows(s_slopes, t_slopes)

# Calculate slopes over time
slopes2 <- slopes %>%
  split(.$id) %>%
  purrr::map(~lm(value ~ year, data = .x)) |>
  purrr::map_df(broom::tidy, .id = 'Year') |>
  filter(term == "year") |> 
  mutate(upr = estimate + 1.96*std.error,
         lwr = estimate - 1.96*std.error,
         sig = ifelse(estimate > lwr & estimate < upr, "sig.", "not sig.")) |> 
  rename(id = Year, 
         year_slope = estimate) |> 
  dplyr::select(id, year_slope, sig, lwr, upr) |> 
  separate(id, into = c("species", "life_stage", "variable", "quarter"), remove = FALSE) |> 
  #mutate(group3 = paste(str_to_title(species), str_to_title(life_stage), paste0("Q", quarter)))
  mutate(group3 = paste(str_to_title(species), str_to_title(life_stage)),
         x = "x") |> 
  mutate(quarter2 = ifelse(quarter == 1, "Quarter 1", "Quarter 4"))
filter: removed 28 rows (50%), 28 rows remaining
mutate: new variable 'upr' (double) with 28 unique values and 0% NA
        new variable 'lwr' (double) with 28 unique values and 0% NA
        new variable 'sig' (character) with one unique value and 0% NA
rename: renamed 2 variables (id, year_slope)
mutate: new variable 'group3' (character) with 7 unique values and 0% NA
        new variable 'x' (character) with one unique value and 0% NA
mutate: new variable 'quarter2' (character) with 2 unique values and 0% NA
hlines <- slopes2 |> 
  filter(species == "Environment")
filter: removed 24 rows (86%), 4 rows remaining
pal <- brewer.pal(name = "Dark2", n = 8)[3:8]

slopes2 |> 
  filter(!species == "Environment") |> 
  ggplot(aes(x, year_slope, ymin = lwr, ymax = upr, color = group3)) + 
  geom_rect(data = hlines, aes(xmin = -Inf, xmax = Inf, ymin = lwr, ymax = upr, fill = "Env. slope 95% CI"),
            inherit.aes = FALSE, alpha = 0.1) +
  geom_hline(yintercept = 0, linetype = 3, color = "tomato2", alpha = 1) +
  geom_hline(data = hlines, aes(yintercept = year_slope, linetype = "Env. slope"), alpha = 0.7) +
  geom_point(size = 2.5, position = position_dodge(width = 0.4)) + 
  geom_errorbar(position = position_dodge(width = 0.4), width = 0) +
  ggh4x::facet_grid2(variable ~ quarter2, scales = "free") +
  scale_fill_manual(values = "grey10") +
  scale_linetype_manual(values = 2) +
  scale_color_manual(values = pal) +
  labs(color = "", y = "Year-slope", linetype = "", fill = "") +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.x = element_blank(), 
        legend.position = "bottom") +
  NULL
filter: removed 4 rows (14%), 24 rows remaining

ggsave(paste0(home, "/figures/year_slopes.pdf"), width = 17, height = 17, units = "cm")

Plot spatially-varying quarter effects, spatial predictions, and ST random fields

grid_preds <- grid_preds |> 
  mutate(group2 = paste(species, life_stage, sep = " ")) 
mutate: new variable 'group2' (character) with 6 unique values and 0% NA
names(grid_preds)
 [1] "X"                 "Y"                 "year"             
 [4] "lon"               "lat"               "depth"            
 [7] "substrate"         "quarter"           "month"            
[10] "month_year"        "oxy"               "temp"             
[13] "sal"               "sub_div"           "ices_rect"        
[16] "temp_sc"           "temp_sq"           "oxy_sc"           
[19] "oxy_sq"            "sal_sc"            "depth_sc"         
[22] "depth_sq"          "quarter_f"         "year_f"           
[25] "est"               "est_non_rf"        "est_rf"           
[28] "omega_s"           "zeta_s_quarter_f1" "zeta_s_quarter_f4"
[31] "epsilon_st"        "group"             "species"          
[34] "life_stage"        "group2"           
# Plot spatially-varying effect
sv <- grid_preds |> 
  #filter(group == pull(grid_preds, group)[1]) |> 
  filter(year == 1999) |> 
  pivot_longer(c("zeta_s_quarter_f1", "zeta_s_quarter_f4"),
               names_to = "field", values_to = "zeta") |> 
  mutate(Quarter = ifelse(field == "zeta_s_quarter_f1", "Quarter 1", "Quarter 4"))
filter: removed 5,465,376 rows (97%), 195,192 rows remaining
pivot_longer: reorganized (zeta_s_quarter_f1, zeta_s_quarter_f4) into (field, zeta) [was 195192x35, now 390384x35]
mutate: new variable 'Quarter' (character) with 2 unique values and 0% NA
plot_map_fc +
  geom_raster(data = sv, aes(X*1000, Y*1000, fill = zeta)) + 
  facet_grid(Quarter ~ group2) +
  scale_fill_gradient2(name = "Zeta")
Warning: Removed 18264 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/sv_effects.pdf"), width = 19, height = 9, units = "cm")
Warning: Removed 18264 rows containing missing values (`geom_raster()`).
for(i in unique(grid_preds$group2)){
  
  dd <- grid_preds |> filter(group2 == i)
  j <- pull(dd, group)[1]
  
  # Quarter 1
  print(
    plot_map_fc +
      geom_raster(data = filter(dd, quarter == 1), aes(X*1000, Y*1000, fill = exp(est))) +
      facet_wrap(~year) + 
      scale_fill_viridis(trans = "sqrt", name = "Biomass density (kg/km)",
                         # trim extreme high values to make spatial variation more visible
                         na.value = "yellow", limits = c(0, exp(quantile(filter(dd, quarter == 4)$est, 0.999)))) +
      labs(caption = paste("maximum estimated biomass density =", round(max(exp(filter(dd, quarter == 4)$est)))),
           title = i, subtitle = "Quarter 1")
    )
    
    ggsave(paste0(home, paste0("/figures/supp/spatial_biomass", "_q1_", j, ".pdf")), width = 17, height = 17, units = "cm")
    
    
  # Quarter 4
  print(
    plot_map_fc +
      geom_raster(data = filter(dd, quarter == 4), aes(X*1000, Y*1000, fill = exp(est))) +
      facet_wrap(~year) + 
      scale_fill_viridis(trans = "sqrt", name = "Biomass density (kg/km)",
                         # trim extreme high values to make spatial variation more visible
                         na.value = "yellow", limits = c(0, exp(quantile(filter(dd, quarter == 4)$est, 0.999)))) +
      labs(caption = paste("maximum estimated biomass density =", round(max(exp(filter(dd, quarter == 4)$est)))),
           title = i, subtitle = "Quarter 4")
    )
    
    ggsave(paste0(home, paste0("/figures/supp/spatial_biomass", "_q4_", j, ".pdf")), width = 17, height = 17, units = "cm")
  
}
filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 4,717,140 rows (83%), 943,428 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining
filter: removed 471,714 rows (50%), 471,714 rows remaining

Warning: Removed 22069 rows containing missing values (`geom_raster()`).
Removed 22069 rows containing missing values (`geom_raster()`).